{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "  import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.sparse as sp\n",
    "from scipy.sparse.linalg import svds\n",
    "from sklearn.metrics.pairwise import pairwise_distances\n",
    "from sklearn.cross_validation import train_test_split\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from math import sqrt\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import warnings; warnings.simplefilter('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "song_df = pd.read_csv('triplet_dataset_sub_song_merged.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>user</th>\n",
       "      <th>song</th>\n",
       "      <th>listen_count</th>\n",
       "      <th>title</th>\n",
       "      <th>release</th>\n",
       "      <th>artist_name</th>\n",
       "      <th>year</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>d6589314c0a9bcbca4fee0c93b14bc402363afea</td>\n",
       "      <td>SOADQPP12A67020C82</td>\n",
       "      <td>12</td>\n",
       "      <td>You And Me Jesus</td>\n",
       "      <td>Tribute To Jake Hess</td>\n",
       "      <td>Jake Hess</td>\n",
       "      <td>2004</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>d6589314c0a9bcbca4fee0c93b14bc402363afea</td>\n",
       "      <td>SOAFTRR12AF72A8D4D</td>\n",
       "      <td>1</td>\n",
       "      <td>Harder Better Faster Stronger</td>\n",
       "      <td>Discovery</td>\n",
       "      <td>Daft Punk</td>\n",
       "      <td>2007</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>d6589314c0a9bcbca4fee0c93b14bc402363afea</td>\n",
       "      <td>SOANQFY12AB0183239</td>\n",
       "      <td>1</td>\n",
       "      <td>Uprising</td>\n",
       "      <td>Uprising</td>\n",
       "      <td>Muse</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>d6589314c0a9bcbca4fee0c93b14bc402363afea</td>\n",
       "      <td>SOAYATB12A6701FD50</td>\n",
       "      <td>1</td>\n",
       "      <td>Breakfast At Tiffany's</td>\n",
       "      <td>Home</td>\n",
       "      <td>Deep Blue Something</td>\n",
       "      <td>1993</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>d6589314c0a9bcbca4fee0c93b14bc402363afea</td>\n",
       "      <td>SOBOAFP12A8C131F36</td>\n",
       "      <td>7</td>\n",
       "      <td>Lucky (Album Version)</td>\n",
       "      <td>We Sing.  We Dance.  We Steal Things.</td>\n",
       "      <td>Jason Mraz &amp; Colbie Caillat</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                       user                song  listen_count  \\\n",
       "0  d6589314c0a9bcbca4fee0c93b14bc402363afea  SOADQPP12A67020C82            12   \n",
       "1  d6589314c0a9bcbca4fee0c93b14bc402363afea  SOAFTRR12AF72A8D4D             1   \n",
       "2  d6589314c0a9bcbca4fee0c93b14bc402363afea  SOANQFY12AB0183239             1   \n",
       "3  d6589314c0a9bcbca4fee0c93b14bc402363afea  SOAYATB12A6701FD50             1   \n",
       "4  d6589314c0a9bcbca4fee0c93b14bc402363afea  SOBOAFP12A8C131F36             7   \n",
       "\n",
       "                           title                                release  \\\n",
       "0               You And Me Jesus                   Tribute To Jake Hess   \n",
       "1  Harder Better Faster Stronger                              Discovery   \n",
       "2                       Uprising                               Uprising   \n",
       "3         Breakfast At Tiffany's                                   Home   \n",
       "4          Lucky (Album Version)  We Sing.  We Dance.  We Steal Things.   \n",
       "\n",
       "                   artist_name  year  \n",
       "0                    Jake Hess  2004  \n",
       "1                    Daft Punk  2007  \n",
       "2                         Muse     0  \n",
       "3          Deep Blue Something  1993  \n",
       "4  Jason Mraz & Colbie Caillat     0  "
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "song_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Train test split\n",
    "n_unique_users = song_df['user'].nunique()\n",
    "n_unique_items = song_df['song'].nunique()\n",
    "\n",
    "train_data, test_data = train_test_split(song_df, test_size=0.4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## User based&Item based"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "users_dictionary = {}\n",
    "items_dictionary = {}\n",
    "ci = 0 \n",
    "cu = 0\n",
    "\n",
    "for i in range(0, song_df.shape[0]-1):\n",
    "    if (items_dictionary.get(song_df.iloc[i].song) == None):\n",
    "        items_dictionary.update({ song_df.iloc[i].song : ci})\n",
    "        ci+=1\n",
    "        \n",
    "    if (users_dictionary.get(song_df.iloc[i].user) == None):\n",
    "        users_dictionary.update({song_df.iloc[i].user : cu})\n",
    "        cu+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "item_val_array = np.fromiter(iter(items_dictionary.values()), dtype=int)\n",
    "item_val_array = np.sort(item_val_array)\n",
    "\n",
    "user_val_array = np.fromiter(iter(users_dictionary.values()), dtype = int)\n",
    "user_val_array = np.sort(user_val_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_data_matrix = np.zeros((n_unique_users, n_unique_items))\n",
    "for line in train_data.itertuples():\n",
    "    train_data_matrix[users_dictionary.get(line[1]), items_dictionary.get(line[2])] = line[3]\n",
    "    \n",
    "test_data_matrix = np.zeros((n_unique_users, n_unique_items))\n",
    "for line in test_data.itertuples():\n",
    "    test_data_matrix[users_dictionary.get(line[1]), items_dictionary.get(line[2])] = line[3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(98632, 800)"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "train_data_matrix.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 相似度计算"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": [
    "user_similarity = pairwise_distances(train_data_matrix, metric = 'cosine')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "item_similarity = pairwise_distances(train_data_matrix.T, metric = 'cosine')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 打分预测"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict_non_biased(ratings, similarity, kind = 'user'):\n",
    "    if kind == 'user':\n",
    "        mean_user_rating = np.nanmean(ratings, axis = 1)\n",
    "        #You use np.newaxis so that mean_user_rating has same format as ratings\n",
    "        ratings_diff = (ratings - mean_user_rating[:, np.newaxis])\n",
    "        pred = mean_user_rating[:, np.newaxis] + similarity.dot(ratings_diff) / np.array([np.abs(similarity).sum(axis = 1)]).T\n",
    "    elif kind == 'item':\n",
    "        pred = ratings.dot(similarity) / np.array([np.abs(similarity).sum(axis = 1)])\n",
    "    return pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "user_prediction_non_biased = predict_non_biased(train_data_matrix, user_similarity, kind = 'user')\n",
    "item_prediction_non_biased = predict_non_biased(train_data_matrix, item_similarity, kind = 'item')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [],
   "source": [
    "#取相似度前20进行预测\n",
    "def predict_topk_non_biased(ratings, similarity, kind='user', k=40):\n",
    "    pred = np.zeros(ratings.shape)\n",
    "    if kind == 'user':\n",
    "        user_bias = ratings.mean(axis=1)\n",
    "        ratings = (ratings - user_bias[:, np.newaxis]).copy()\n",
    "        for i in range(ratings.shape[0]):\n",
    "            top_k_users = [np.argsort(similarity[:,i])[:-k-1:-1]]\n",
    "            for j in range(ratings.shape[1]):\n",
    "                pred[i, j] = similarity[i, :][top_k_users].dot(ratings[:, j][top_k_users]) \n",
    "                pred[i, j] /= np.sum(np.abs(similarity[i, :][top_k_users]))\n",
    "        pred += user_bias[:, np.newaxis]\n",
    "    if kind == 'item':\n",
    "        item_bias = ratings.mean(axis=0)\n",
    "        ratings = (ratings - item_bias[np.newaxis, :]).copy()\n",
    "        for j in range(ratings.shape[1]):\n",
    "            top_k_items = [np.argsort(similarity[:,j])[:-k-1:-1]]\n",
    "            for i in range(ratings.shape[0]):\n",
    "                pred[i, j] = similarity[j, :][top_k_items].dot(ratings[i, :][top_k_items].T) \n",
    "                pred[i, j] /= np.sum(np.abs(similarity[j, :][top_k_items])) \n",
    "        pred += item_bias[np.newaxis, :]\n",
    "        \n",
    "    return pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [],
   "source": [
    "topk_user_pred_non_biased = predict_topk_non_biased(train_data_matrix, user_similarity, kind = 'user', k = 20)\n",
    "topk_item_pred_non_biased = predict_topk_non_biased(train_data_matrix, item_similarity, kind = 'item', k = 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mse(pred, actual):\n",
    "    # Ignore nonzero terms.\n",
    "    pred = pred[actual.nonzero()].flatten()\n",
    "    actual = actual[actual.nonzero()].flatten()\n",
    "    return mean_squared_error(pred, actual)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "User-based CF MSE: 186.3150485566105\n",
      "Item-based CF MSE: 187.63338250929647\n"
     ]
    }
   ],
   "source": [
    "print('User-based CF MSE: ' + str(mse(user_prediction_non_biased, test_data_matrix)))\n",
    "print('Item-based CF MSE: ' + str(mse(item_prediction_non_biased, test_data_matrix)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Top K User-based CF MSE: ' + str(mse(topk_user_pred_non_biased, test_data_matrix)))\n",
    "print('Top K Item-based CF MSE: ' + str(mse(topk_item_pred_non_biased, test_data_matrix)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SVD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [],
   "source": [
    "#adjusting parameter K for SVD\n",
    "error_rate = []\n",
    "\n",
    "# Will take some time\n",
    "for i in range(1,50):  \n",
    "    u, s, vt = svds(train_data_matrix, k = i)\n",
    "    s_diag_matrix=np.diag(s)\n",
    "    X_pred = np.dot(np.dot(u, s_diag_matrix), vt)\n",
    "    error_rate.append(mse(X_pred, test_data_matrix))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(50,70):  \n",
    "    u, s, vt = svds(train_data_matrix, k = i)\n",
    "    s_diag_matrix=np.diag(s)\n",
    "    X_pred = np.dot(np.dot(u, s_diag_matrix), vt)\n",
    "    error_rate.append(mse(X_pred, test_data_matrix))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA34AAAJQCAYAAADR+LbmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmcXWV9P/DPCSEsmQAqQQg7JA0CQZSIW91wQ0EEXKKirfwqWMBalVaLddeqdV9RW0HrT+HnAlFwQ6yCVqQVBBUlAaIQIATCIkv25fz+eGaaELJMknvuuXPn/X69zuvcOXPnzjejzNzPeZ7n+1R1XQcAAID+NabtAgAAAGiW4AcAANDnBD8AAIA+J/gBAAD0OcEPAACgzwl+AAAAfU7wAwAA6HOCHwAAQJ8T/AAAAPrc2LYL2BI777xzvc8++7RdBgAAQCuuvPLKO+u6nrix543o4LfPPvvkiiuuaLsMAACAVlRVddNwnmeqJwAAQJ8T/AAAAPqc4AcAANDnBD8AAIA+J/gBAAD0OcEPAACgzwl+AAAAfU7wAwAA6HOCHwAAQJ8T/AAAAPqc4AcAANDnBD8AAIA+J/gBAAD0OcEPAACgzwl+AAAAfU7wAwAA6HOCHwAAQJ8T/AAAAPqc4AcAANDnBD8AAIA+J/gBAAD0OcEPAACgzwl+nfapTyUHHpjUdduVAAAAJBH8mnHttcmCBW1XAQAAkETw67zJk8v5hhvarQMAAGCQ4Ndpgh8AANBjBL9O22efZMwYwQ8AAOgZgl+njRtXwp/gBwAA9AjBrwmTJyfXX992FQAAAEkEv2YMBT9bOgAAAD1A8GvC5MnJvfcmd9/ddiUAAACCXyN09gQAAHqI4NeEKVPKWfADAAB6gODXhH33TapKgxcAAKAnCH5N2GabZK+9jPgBAAA9QfBryuTJgh8AANATBL+mCH4AAECPEPyaMnlyctddyT33tF0JAAAwygl+TdHZEwAA6BGCX1Ps5QcAAPQIwa8p++1XzoIfAADQMsGvKdttl+yxh+AHAAC0TvBrks6eAABADxD8mjRlSnL99W1XAQAAjHKCX5MmT04WLEjuvbftSgAAgFFM8GvSUGfPOXParQMAABjVBL8m2dIBAADoAYJfk/bfv5wFPwAAoEWCX5PGj08mTRL8AACAVgl+TZs8WWdPAACgVYJf0+zlBwAAtEzwa9rkycn8+ckDD7RdCQAAMEoJfk2zpQMAANAywa9ptnQAAABaJvg1bSj4afACAAC0RPBr2oQJySMfacQPAABojeDXDTp7AgAALRL8ukHwAwAAWiT4dcPkycmttyaLFrVdCQAAMAoJft0wZUo529IBAABogeDXDbZ0AAAAWiT4dcP++5ez4AcAALRA8OuGnXZKdt5Z8AMAAFoh+HWLzp4AAEBLBL9umTw5uf76tqsAAABGIcGvW6ZMSW6+OVm8uO1KAACAUUbw65ahzp5/+lO7dQAAAKOO4NcttnQAAABaIvh1i+AHAAC0RPDrloc/PHnYwzR4AQAAuk7w66YpU4z4AQAAXSf4dZO9/AAAgBYIft00eXIyd26ydGnblQAAAKOI4NdNkycnq1YlN97YdiUAAMAoIvh101BnTw1eAACALhL8usmWDgAAQAsEv27aeedkxx0FPwAAoKsEv26qKp09AQCArhP8uk3wAwAAuqyx4FdV1dlVVd1RVdU1a1w7tKqqy6uqurqqqiuqqjp88PrDqqqaWVXVb6uq+p+qqg5uqq7WTZ5cunouX952JQAAwCjR5Ijfl5Mcuda1DyV5d13XhyZ5x+DHSfLWJFfXdX1Ikr9K8skG62rX5MnJypW2dAAAALqmseBX1/XPkty99uUkOww+3jHJvMHHByb5z8Gvm5Vkn6qqHtlUba2aMqWcTfcEAAC6ZGyXv98bklxUVdVHUkLnkwav/ybJ8Un+a3D6595J9khye5fra54tHQAAgC7rdnOXU5K8sa7rPZO8MclZg9c/mORhVVVdneTvklyVZMW6XqCqqpMH1wdesWDBgm7U3Fm77JIMDAh+AABA13Q7+P11kvMHH38zyeFJUtf1fXVdnzi49u+vkkxM8qd1vUBd1/9W1/X0uq6nT5w4sRs1d5YtHQAAgC7rdvCbl+Rpg4+PSHJ9klRVtVNVVeMGr78myc/qur6vy7V1z+TJyfXXt10FAAAwSjS2xq+qqnOTPD3JzlVV3ZLknUlOSvLJqqrGJlmS5OTBpz8qyVeqqlqZ5A9J/qapunrC5MnJt7+drFiRjO32MksAAGC0aSx11HX98vV86rB1PPeXSaY0VUvPmTKlhL65c5P99mu7GgAAoM91e6onic6eAABAVwl+bRD8AACALhL82rDbbsl222nwAgAAdIXg1wZbOgAAAF0k+LVlyhTBDwAA6ArBry2TJyd//GOycmXblQAAAH1O8GvL5MnJsmXJLbe0XQkAANDnBL+2DHX21OAFAABomODXFls6AAAAXSL4tWX33ZNttxX8AACAxgl+bRkzJtl/f8EPAABonODXJnv5AQAAXSD4tWny5GTOnGTVqrYrAQAA+pjg16bJk5MlS2zpAAAANErwa9Puu5fz7be3WwcAANDXBL82TZhQzgsXtlsHAADQ1wS/Ng0MlPMDD7RbBwAA0NcEvzYJfgAAQBcIfm0aP76cBT8AAKBBgl+bjPgBAABdIPi1yYgfAADQBYJfm8aNK4fgBwAANEjwa9vAgOAHAAA0SvBr28CAffwAAIBGCX5tM+IHAAA0TPBrm+AHAAA0TPBrm+AHAAA0TPBrm+AHAAA0TPBrm+AHAAA0TPBrm+AHAAA0TPBrm+AHAAA0TPBr29A+fqtWtV0JAADQpwS/tg0MJHWdLF7cdiUAAECfEvzaNjBQzqZ7AgAADRH82jZ+fDkLfgAAQEMEv7YZ8QMAABom+LVN8AMAABom+LVN8AMAABom+LVN8AMAABom+LVtKPgtXNhuHQAAQN8S/NpmxA8AAGiY4Nc2wQ8AAGiY4Ne27bZLqkrwAwAAGiP4ta2qyqif4AcAADRE8OsFgh8AANAgwa8XCH4AAECDBL9eIPgBAAANEvx6geAHAAA0SPDrBYIfAADQIMGvF4wfL/gBAACNEfx6gRE/AACgQYJfLxD8AACABgl+vUDwAwAAGiT49YKBgWT58mTZsrYrAQAA+pDg1wsGBsp54cJ26wAAAPqS4NcLhoKf6Z4AAEADBL9eIPgBAAANEvx6geAHAAA0SPDrBYIfAADQIMGvFwh+AABAgwS/XiD4AQAADRL8eoHgBwAANEjw6wWCHwAA0CDBrxeMH1/Ogh8AANAAwa8XjB2bbLut4AcAADRC8OsV48cLfgAAQCMEv14xMCD4AQAAjRD8eoXgBwAANETw6xWCHwAA0BDBr1cMDCQLF7ZdBQAA0IcEv15hxA8AAGiI4NcrBD8AAKAhgl+vEPwAAICGCH69QvADAAAaIvj1iqHmLqtWtV0JAADQZwS/XjEwUM6LFrVbBwAA0HcEv14xFPxM9wQAADpM8OsVgh8AANAQwa9XCH4AAEBDBL9eIfgBAAANEfx6xfjx5Sz4AQAAHSb49QojfgAAQEMEv14h+AEAAA0R/HqF4AcAADRE8OsVgh8AANAQwa9XbLttMmZMsnBh25UAAAB9RvDrFVVVRv2M+AEAAB0m+PUSwQ8AAGiA4NdLBD8AAKABgl8vEfwAAIAGCH69RPADAAAaIPj1EsEPAABogODXSwQ/AACgAYJfLxH8AACABgh+vUTwAwAAGiD49ZLx40vwq+u2KwEAAPqI4NdLBgaSFSuSZcvargQAAOgjgl8vGRgoZ9M9AQCADhL8eongBwAANEDw6yWCHwAA0ADBr5cMBb+FC9utAwAA6CuCXy8x4gcAADRA8Oslgh8AANAAwa+XCH4AAEADGgt+VVWdXVXVHVVVXbPGtUOrqrq8qqqrq6q6oqqqwwev71hV1YVVVf2mqqrfV1V1YlN19TTBDwAAaECTI35fTnLkWtc+lOTddV0fmuQdgx8nyWlJ/lDX9aOTPD3JR6uqGtdgbb1J8AMAABrQWPCr6/pnSe5e+3KSHQYf75hk3hrXJ1RVVSUZGPy6FU3V1rPGjy9nwQ8AAOigsV3+fm9IclFVVR9JCZ1PGrz+mSQXpATBCUlm1HW9al0vUFXVyUlOTpK99tqr8YK7aqutku22E/wAAICO6nZzl1OSvLGu6z2TvDHJWYPXn5vk6iSTkhya5DNVVe2wrheo6/rf6rqeXtf19IkTJ3aj5u4aGBD8AACAjup28PvrJOcPPv5mksMHH5+Y5Py6uCHJn5Ic0OXaeoPgBwAAdFi3g9+8JE8bfHxEkusHH89N8swkqarqkUmmJvljl2vrDYIfAADQYY2t8auq6tyUDp07V1V1S5J3JjkpySerqhqbZEkG1+oleW+SL1dV9bskVZK31HV9Z1O19bTx4wU/AACgoxoLfnVdv3w9nzpsHc+dl+Q5TdUyohjxAwAAOqzbUz3ZGMEPAADoMMGv1wh+AABAhwl+vUbwAwAAOkzw6zUDA8nChW1XAQAA9BHBr9cMDCSLFiUrV7ZdCQAA0CcEv14zMFDOixa1WwcAANA3BL9eMxT8rPMDAAA6RPDrNYIfAADQYYJfrxH8AACADhP8eo3gBwAAdJjg12sEPwAAoMMEv14j+AEAAB0m+PUawQ8AAOgwwa/XjB9fzoIfAADQIYJfrzHiBwAAdJjg12u22SbZaivBDwAA6BjBr9dUVRn1E/wAAIAOEfx6keAHAAB0kODXiwYGkoUL264CAADoE4JfLzLiBwAAdJDg14sEPwAAoIMEv14k+AEAAB0k+PUiwQ8AAOggwa8XCX4AAEAHCX69SPADAAA6SPDrRUPBr67brgQAAOgDgl8vGhhIVq5Mli5tuxIAAKAPCH69aGCgnE33BAAAOkDw60WCHwAA0EGCXy8aP76cBT8AAKADBL9eZMQPAADoIMGvFwl+AABABwl+vUjwAwAAOkjw60VDwW/hwnbrAAAA+oLg14uM+AEAAB0k+PUiwQ8AAOggwa8Xbb99OQt+AABABwh+vWirrUr4E/wAAIAOEPx61cCA4AcAAHSE4NerBD8AAKBDBL9eJfgBAAAdIvj1KsEPAADoEMGvVwl+AABAhwh+vUrwAwAAOkTw61Xjxwt+AABARwh+vcqIHwAA0CHDCn5VVf1lVVUnDj6eWFXVvs2WheAHAAB0ykaDX1VV70zyliRnDF7aOslXmyyKlOC3ZEmyYkXblQAAACPccEb8jktyTJKFSVLX9bwkE5osipTglyQLF7ZbBwAAMOINJ/gtq+u6TlInSVVV45stiSSCHwAA0DHDCX7fqKrqC0l2qqrqpCQ/TvLFZsvif4OfdX4AAMAWGruxJ9R1/ZGqqp6d5L4kU5O8o67rixuvbLQT/AAAgA7ZaPCrqupf67p+S5KL13GNpgh+AABAhwxnquez13HteZ0uhLUIfgAAQIesd8SvqqpTkpyaZL+qqn67xqcmJPlF04WNeoIfAADQIRua6nlOkh8k+UCSf1rj+v11Xd/daFUIfgAAQMesN/jVdX1vknuTvDxJqqraJcm2SQaqqhqo63pud0ocpQQ/AACgQza6xq+qqhdUVXV9kj8luTTJjSkjgTRJ8AMAADpkOM1d3pfkCUmuq+t63yTPjDV+zRs3Ltl6a8EPAADYYsMJfsvrur4ryZiqqsbUdf3TJIc2XBdJMn684AcAAGyxje7jl+TPVVUNJPlZkq9VVXVHkhXNlkWSMt1T8AMAALbQcEb8XphkUZI3JvlhkjlJXtBkUQwS/AAAgA7Y4IhfVVVbJflOXdfPSrIqyX90pSoKwQ8AAOiADY741XW9Msmiqqp27FI9rEnwAwAAOmA4a/yWJPldVVUXJ1k4dLGu69c3VhXFwEAyb17bVQAAACPccILf9wYPus2IHwAA0AEbDX51XVvX15Ymgt+iRcncuauPm24q56c+Nfmbv+ns9wIAAHrCcEb8aMuWBL/LL08uu+zB4W7u3OTOOx/8vDFjkgkTknPOSZ74xOTAA7e8bgAAoKcIfr1sKPjVdVJVw/+6VauS5zwnuf/+8hp7753stVfyuMetfrzXXuXxpEnJPfckU6cmp56a/PSnm/a9AACAnjec7Rw+WNf1P3apHtY0MFBC3JIlyXbbDf/rbrmlhL5Pfzo57bSNB7mJE5MPfCD5278tI38nnLBldQMAAD1lONs5HFZVhoBaMTBQzps63XP27HKeNm34o3eveU1y+OHJ6acn9967ad8PAADoaRsMfoOuSvKdqqpeVVXV8UNH04WRzQ9+s2aV89Spw/+arbZKzjwzueOO5O1v37TvBwAA9LThBL+HJ7kryRFJXjB4HN1kUQzakhG/HXZIHvnITfu6ww5LTjkl+exnk6uu2rSvBQAAetZwtnM4sRuFsA5bMuJ3wAGb16Tlfe9LvvWtEgAvu6x0/QQAAEa0jb6rr6pqj6qqZlZVdUdVVbdXVXVeVVV7dKO4UW9LRvw2ZZrnmh72sOTDH07++7+Ts8/evNcAAAB6ynCGc76U5IIkk5LsnuTCwWs0bfz4ct6U4PfAA6Wr5wEHbP73fdWrkqc8JXnLWx667x8AADDiDCf4Tazr+kt1Xa8YPL6cZGLDdZFs3ojfddeV8+aO+CVliuiZZ5bunmecsfmvAwAA9IThBL87q6p6ZVVVWw0er0xp9kLTNif4DW3lsCXBL0kOPjh5wxuSL34xufzyLXstAACgVcMJfv8nyUuTzE9yW5IXD16jaZsT/GbNKg1ZJk/e8u//zncmu+9eGr2sWLHlrwcAALRig8Gvqqqtkryorutj6rqeWNf1LnVdH1vX9U1dqm902377cl64cPhfM3t2ss8+ybbbbvn3nzAh+fjHk6uvTj73uS1/PQAAoBUbDH51Xa9M8sIu1cLaxowpDV42dcRvSxq7rO3FL06e/ezkbW9Lbrutc68LAAB0zXCmev6iqqrPVFX1lKqqHjt0NF4ZxcDA8IPfqlWlucuWru9bU1Uln/lMsmRJ8o//2LnXBQAAumajG7gnedLg+T1rXKuTHNH5cniITQl+t9ySLF7c2eCXJH/xF8mb31w2d3/Na5KnP72zrw8AADRqg8GvqqoxST5X1/U3ulQPa9uU4DdrVjl3cqrnkDPOSL761eTUU8uav3HjOv89AACARmxsjd+qJK/rUi2sy6YEv05t5bAu22+ffPrTybXXJp/4ROdfHwAAaMxw1vhdXFXVP1RVtWdVVQ8fOhqvjGJTR/x23DF55CObqeXoo5Njjkne9a7Vo4sAAEDPG+4+fqcl+VmSKwePK5osijVs6ojf1KmlIUtTPve5Mvr3ilcky5Y1930AAICO2Wjwq+t633Uc+3WjOLLpI35NrO9b06RJyVlnJVddVbZ4AAAAet56g19VVW9e4/FL1vrc+5ssijUMN/g98EBy663NrO9b2wtfmLz2tcmHP5z85382//0AAIAtsqERv5et8fiMtT53ZAO1sC7DDX7XXVfO3Qh+SfKxj5Xv9Vd/ldx1V3e+JwAAsFk2FPyq9Txe18c0Zfz4ZOnSZPnyDT+vya0c1mX77ZNzz00WLEhOOimp6+58XwAAYJNtKPjV63m8ro9pysBAOS9cuOHnzZ6djBmTTJ7cfE1DHvOY5P3vT2bOTL74xe59XwAAYJNsKPg9uqqq+6qquj/JIYOPhz6e1qX6GAp+G5vuOWtWsu++yTbbNF/Tmt70puSZz0ze8IbV+wgCAAA9Zb3Br67rreq63qGu6wl1XY8dfDz08dbdLHJUG27wG9rKodvGjEn+4z+SbbdNTjjBFg8AANCDhrOPH20aTvBbtao0d2kj+CXJ7ruXLR6uvDJ5+9vbqQEAAFgvwa/XDWeN3803J4sXd6+xy7oce2xy8slli4ef/GTTv76ukzlzNIkBAIAGCH69bjgjfkNr69oa8RvysY8lf/EXm7bFw333JWeemUybVhrTnLH2ziEAAMCWEvx63XCCX7e3clif8eOTc85J7rijjP5taPTummuS004r00RPO600pXnBC5J//dfk3/+9ezUDAMAoIPj1uuGO+O24Y7LLLt2paUMe+9jkX/4lOf/85OyzH/y5ZcuSb3wjedrTygjfWWclxx+f/Pd/J1dcUb7myCOTU05JLr64nfoBAKAPCX69brgjflOnJlXVnZo25vTTkyOOSF7/+tJ05pZbkne8I9l772TGjLIm8UMfKtf/4z+Sww8vtY8dm3z968mBByYvfnEZFQQAALbY2LYLYCPGjy/njY34PfOZ3alnOMaMSb7yleSQQ5KnPKWs91u1Knn+88u0zuc+tzxnXXbYIfne95LHPz456qgyGrjrrt2tHwAA+owRv143blw51hf87r8/ufXW9hu7rG333cto3oQJZQRwzpzku99Nnve89Ye+IXvumVx4YXLnnWXd36JF3akZAAD6lOA3EgwMrD/4XXddObfd2GVdjj46ueGG0rBl33037WsPOyw599yyN+ArX5msXNlMjQAAMAo0Fvyqqjq7qqo7qqq6Zo1rh1ZVdXlVVVdXVXVFVVWHD17/x8FrV1dVdU1VVSurqnp4U7WNOBsKfr2ylUMTjjkm+cQnkpkzkze/ue1qAABgxGpyxO/LSY5c69qHkry7rutDk7xj8OPUdf3huq4PHbx+RpJL67q+u8HaRpYNBb9Zs8rUycmTu1tTt7z+9cnf/V3ZI/DMM9uuBgAARqTGmrvUdf2zqqr2Wftykh0GH++YZN46vvTlSc5tqq4RaWMjfvvuW/bB61cf/3jypz+VALjPPqVJDAAAMGzdXuP3hiQfrqrq5iQfSRnd+19VVW2fMkp4Xpfr6m3jx294xK8fp3muaautynq/Rz+6bAfxm9+0XREAAIwo3Q5+pyR5Y13XeyZ5Y5Kz1vr8C5L8YkPTPKuqOnlwfeAVCxYsaLDUHrK+Eb9Vq5Lrr+/Nxi6dNjBQOn3uuGPZ5uHWW9uuCAAARoxuB7+/TnL+4ONvJjl8rc+/LBuZ5lnX9b/VdT29ruvpEydObKDEHrS+4Hfzzcnixf0/4jdk993LHn/33lu2ebjyyrLJ++zZyR//WH4e8+eXfQPvu6/8bFasSOq67coBAKBV3d7AfV6SpyW5JMkRSa4f+kRVVTsOfu6VXa6p960v+M2aVc6jYcRvyKMfnXzjG2WriOnTh/91b3pT8tGPNlcXAAD0sMaCX1VV5yZ5epKdq6q6Jck7k5yU5JNVVY1NsiTJyWt8yXFJflTX9cKmahqxBgaShev4sfTzVg4b8rznlXV+11+fLF9eRvWWL1//8c1vlpFCwQ8AgFGqya6eL1/Ppw5bz/O/nLIFBGsbGvGr66SqVl+fNausedtll/Zqa8vBB5djOFasSN73vjL1c7vtmq0LAAB6ULfX+LE5BgZK6Fu8+MHXZ88u0zzXDIM81LRppRHOtde2XQkAALRC8BsJBgbKee11fqNhK4dOmDatnH/3u3brAACAlgh+I8G6gt/99yfz5o2uxi6ba/LkZNttBT8AAEYtwW8kWFfwu+66cjbit3FbbZUceKDgBwDAqCX4jQTrCn5DWzkIfsMzbZrgBwDAqCX4jQTrCn6zZydjxpRpjGzctGnJbbeVzd0BAGCUEfxGgvWN+O27b7LNNu3UNNJo8AIAwCgm+I0E6xvx09hl+IaC329/224dAADQAsFvJFg7+K1aVZq7WN83fLvumjziEUb8AAAYlQS/kWD8+HIeCn5z5yZLlgh+m6KqNHgBAGDUEvxGgu22K8FlKPjNnl3OpnpummnTkmuuKSOmAAAwigh+I8GYMWXUbyj42cph80yblixcmNx4Y9uVAABAVwl+I8XAwINH/HbaKdlll3ZrGml09gQAYJQS/EaKgYEyWpWUEb+pU8v0T4bvoIPKWfADAGCUEfxGirVH/Kzv23QTJpS9DwU/AABGGcFvpBgKfvffn8ybZ33f5jrkEMEPAIBRR/AbKYaC31BHT8Fv80ybVvZAXLq07UoAAKBrBL+RYu3gZ6rn5pk2LVm5Mrn22rYrAQCArhH8Roqh4DdrVtneYf/9265oZNLZEwCAUUjwGynWHPHbb79km23armhkmjKl/OwEPwAARhHBb6RYc8TP+r7NN3Zs8qhHCX4AAIwqgt9IMTCQLFtWRvwEvy0zbZrgBwDAqCL4jRQDA+W8bJnGLltq2rTk1luTe+5puxIAAOgKwW+kGAp+iRG/LaXBCwAAo4zgN1KsGfyM+G0ZwQ8AgFFG8Bspxo8v5512SiZObLeWkW7SpORhD0t++9u2KwEAgK4Q/EaKoRG/qVOTqmq3lpGuqjR4AQBgVBH8Roqh4GeaZ2dMm5Zcc01S121XAgAAjRP8Roo1R/zYctOmJfffn9x0U9uVAABA4wS/kWLffZPjjkuOPbbtSvqDBi8AAIwigt9Isc02yfnnJ496VNuV9IeDDy5nwQ8AgFFA8GN02mGHZO+9BT8AAEYFwY/RS2dPAABGCcGP0WvatGT27GTZsrYrAQCARgl+jF7TpiUrViSzZrVdCQAANErwY/TS2RMAgFFC8GP0mjo12XprwQ8AgL4n+DF6bb11csABgh8AAH1P8GN009kTAIBRQPBjdJs2Lbn55uTPf267EgAAaIzgx+h2yCHlfM017dYBAAANEvwY3XT2BABgFBD8GN322CPZcUfBDwCAvib4MbpVlQYvAAD0PcEPhoJfXbddCQAANELwg2nTknvvLd09AQCgDwl+oMELAAB9TvCDgw8uZ8EPAIA+JfjBTjsle+4p+AEA0LcEP0h09gQAoK8JfpCU4DdrVrJ8eduVAABAxwl+kJTgt3x5Mnt225UAAEDHCX6Q6OwJAEBfE/wgSQ44IBk7VvADAKAvCX6QJOPGJVOnCn4AAPQlwQ+G6OwJAECfEvxgyLRpyU03Jffd13YlAADQUYIfDBlq8HLNNe3WAQAAHSb4wRCdPQEA6FOCHwzZe+9kwgTBDwCAviP4wZCqSg4+WPADAKDvCH6wpic+MfnFL5If/KDtSgAAoGMEP1jTu95V1vq95CXJr3/ddjUAANARgh+sacKE5HvfSx7+8OSoo8r2DgAAMMIlkMc9AAAgAElEQVQJfrC2SZPKVM/Fi5PnPz/585/brggAALaI4AfrctBBycyZyfXXJ8cdlyxd2nZFAACw2QQ/WJ9nPCP50peSSy5J/s//Seq67YoAAGCzjG27AOhpJ5xQ1vn98z8n++yT/Mu/tF0RAABsMsEPNuaMM5Ibb0ze//6yyfvJJ7ddEQAAbBLBDzamqpIzz0xuvTU55ZRkjz1K0xcAABghrPGD4Rg7Nvn615NDD01e+tLkyivbrggAAIZN8IPhGhhIvvvdZOedk6OPLtM/AQBgBBD8YFPstlvy/e8nS5aU6Z733NN2RQAAsFGCH2yqAw9Mvv3tZM6cssffypVtVwQAABsk+MHmeNrTkk9+Mrn00uTnP2+7GgAA2CDBDzbXq16VbL99afoCAAA9TPCDzTV+fPKCFyTnnZesWNF2NQAAsF6CH2yJGTOSBQuSn/607UoAAGC9BD/YEs97XjJhgumeAAD0NMEPtsS22yYvfGFy/vnJsmVtVwMAAOsk+MGWmjGj7Of34x+3XQkAAKyT4Adb6jnPSXbayXRPAAB6luAHW2rcuLKR+7e/nSxZ0nY1AADwEIIfdMKMGcl99yUXXdR2JQAA8BCCH3TCEUckj3iE6Z4AAPQkwQ86Yeutkxe9KLnggmTRorarAQCABxH8oFNmzEgWLky+//22KwEAgAcR/KBTnvrUZJddTPcEAKDnCH7QKWPHJi9+cfK97yUPPNB2NQAA8L8EP+ikGTOSxYuTCy9suxIAAPhfgh900l/+ZTJpkumeAAD0FMEPOmnMmOQlL0l+8IPk3nvbrgYAAJIIftB5M2Yky5Yl3/lO25UAAEASwQ867wlPSPbay3RPAAB6huAHnVZVyUtfmvzoR8ndd7ddDQAACH7QiBkzkhUrkpkz264EAAAEP2jEYYcl++1nuicAAD1B8IMmVFUZ9fvJT5IFC9quBgCAUU7wg6bMmJGsXJmcf37blQAAMMoJftCUQw5Jpk413RMAgNYJftCUoemel16azJ/fdjUAAIxigh80acaMZNWq5FvfarsSAABGMcEPmnTggcnBB5vuCQBAqxoLflVVnV1V1R1VVV2zxrVDq6q6vKqqq6uquqKqqsPX+NzTB6//vqqqS5uqC7puxozkv/4rueWWtisBAGCUanLE78tJjlzr2oeSvLuu60OTvGPw41RVtVOSM5McU9f1QUle0mBd0F0zZpTzN7/Zbh0AAIxajQW/uq5/luTutS8n2WHw8Y5J5g0+fkWS8+u6njv4tXc0VRd03ZQpyWMeY7onAACtGdvl7/eGJBdVVfWRlND5pMHrf5Fk66qqLkkyIckn67r+Spdrg+bMmJH80z8lb3pTWfc3dWo5Jk4s3T8BAKBB3Q5+pyR5Y13X51VV9dIkZyV51mAdhyV5ZpLtkvyyqqrL67q+bu0XqKrq5CQnJ8lee+3VtcJhi7zqVWWq55lnJkuXrr7+sIetDoFTpyYHHFDO+++fbLNNe/UCANBXqrqum3vxqtonyXfruj548ON7k+xU13VdVVWV5N66rneoquqfkmxb1/W7Bp93VpIf1nW9wUVR06dPr6+44orG6oeOW7kymTs3mT07mTWrnIeOefNWP2/s2OSLX0z++q/bqxUAgJ5XVdWVdV1P39jzuj3iNy/J05JckuSIJNcPXv9Oks9UVTU2ybgkj0/y8S7XBs3baqtk333LceRavY/uuy+57roSAj/zmeTv/z557nOTXXdtp1YAAPpGY8Gvqqpzkzw9yc5VVd2S5J1JTkryycGAtySDUzbrur62qqofJvltklVJvljX9TXrfGHoVzvskEyfXo7HPS6ZNi35h39IvvrVtisDAGCEa3SqZ9NM9aSvvf3tyfvel/zkJ8kzntF2NQAA9KDhTvVsch8/YEu89a1lSuippybLlrVdDQAAI5jgB71qu+3KWr9Zs5KPfrTtagAAGMEEP+hlz39+cvzxyXvfm9x4Y9vVAAAwQgl+0Os+8YlkzJjk9a9vuxIAAEYowQ963Z57Ju98Z3LhhckFF7RdDQAAI5DgByPBG96QHHRQ8nd/lyxc2HY1AACMMIIfjARbb5187nPJ3LlliwcAANgEgh+MFE95SvLqVycf+Ujyhz+0XQ0AACOI4AcjyYc+lEyYkJx2WlLXbVcDAMAIIfjBSDJxYvLBDyaXXJJ87WttVwMAwAgh+MFI85rXJIcfnpx+enLPPW1XAwDACCD4wUgzZkzy+c8nd96ZvO1tbVcDAMAIIPjBSPSYxySve13p9PmrX7VdDQAAPU7wg5HqPe9JHvnI5JRTkpUr264GAIAeJvjBSLXjjsnHP55ceWUZ+QMAgPUQ/GAkmzEjec5zkre8JbnhhrarAQCgRwl+MJJVVXLWWcm4ccmrXpWsWNF2RQAA9CDBD0a6PfZIzjwzufzy5F//te1qAADoQYIf9IOXvzx52cuSd72rrPkDAIA1CH7QLz772dLl85WvTBYvbrua0U2XVQCgxwh+0C8e/vDkS19KZs1Kzjij7WpGr5tvTnbaKTn33LYrAQD4X4If9JNnPzv5u79LPvnJ5Mc/brua0em885IHHkj+4R+ShQvbrgYAIIngB/3ngx9Mpk5NXv3q5J572q5m9Jk5M9l552TevOSjH227GgCAJIIf9J/tt0+++tXk9tuT172u7WpGlwULkv/6r+SUU5IXv7h0WZ03r+2qAAAEP+hL06cn73hHcs45yde/vnmvsWhRMnduZ+vqdxdemKxalRx3XBl5XbEiedvb2q4KAEDwg751xhnJ4x9fRp9uvXX4X7dkSVkjuN9+ZcroggXN1dhvZs5M9t47OfTQZP/9y3rLL385ufrqtisDAEY5wQ/61dixyVe+kixdmpx4YhmJ2pBly5LPfS6ZPDl5wxvKxvBLliQXXNCdeke6++9PLr64jPZVVbn2treVbqunn57Udbv1Ad0zf35y333t1nDPPcl3v5vccUe7dQA9Q/CDfvYXf1EajFx8cXLmmet+zvLlyVlnleeeemqyzz7JT36S/OpX5fHMmd2seOS66KISso89dvW1nXZK3vnO8vP83vfaqw1ozm23lYD1rnclL3hBMmlSsttuya67Jied1N0R/8WLk29+s9yA2nXX1fU873nJ175WOg4Do1ZVj+C70NOnT6+vuOKKtsuA3lbXyVFHJZdckvz618kBB5TrK1eWNYDvfncyZ07yuMcl731v8pznrB6xOv305DOfKdM9d9ihtX/CiHDCCcmPflTu9G+11erry5cnBx9cfqa/+12y9dbt1QhsmfnzkyuvTK64opyvvHJ1A6eqSh71qOSww5LHPjb5wx9Ko63Fi5MnPzk57bTkRS9Kxo3rbE0rVpSbS+eck5x/fpl9sNtuycteljz3ueV3/znnlDXb229fbk6dcELZ/sfvI+gLVVVdWdf19I0+T/CDUeC220r42G+/0nXy298uI1GzZ5f1aO95T3L00asD35Bf/CL5y78sm5G/7GXt1D4SLFuW7LJLeVN31lkP/fyFFybHHFNC9Gmndb8+GA0WLiwj7//1XyUMrWl973XquozUL1lSjsWLVz9e+9qiRcmf/1y+rqrKTbTDDivH9Onld+nAwINf/557ki99qcy4mDMneeQjk5NPTl772mT33Tf/31rXyf/8z+oGXrffnuy4Y/kd9IpXJE9/+oNvQK1aVX4uX/taGRG8555k4sRkxowSAh//+If+/gdGDMEPeLDzzitbDEycWEbwDjqojPYdd1wyZj2zvletKtOEnvrU5Bvf6G69I8lFFyVHHlkC3tFHP/TzdZ0885nJb3+b3HBDmQIKbLm77y7TLGfOLP8dLl6cbLddsu22637+usLNttuu/xh6raFj//1Xh7wJE4Zf56pVpb7Pfjb5/vfL79zjjitb7jz1qesPXStWJHfeWX5n33FHOf/+9+Vm3Jw5yTbblN85J5xQpnOu79+9pqVLkx/+sITACy8soXb//UtgPPbY8m9b398EoCcJfsBDnXZa8vOfJ299a/KSlzz4jvD6/O3flulKCxaUN0E81CmnJP/3/5Y3aOt743XVVWVk4PTTkw9/uLv1QT+59dYya2HmzDKNceXK0ozq2GNLmHrqU0tzq141Z05ppHX22WXk7eCDk+c/v4wmLliw+rjjjvL5tY0ZkxxxRAl7xx1XRvo21333lemhX/1qmS5a18nOO5cbVc9+djn22mvzXx/oCsEP6Iwf/aisE7nggtIogAdbtaq86Xzyk8sUqg058cQyNevaa8u0W2B4rruuBL3zzy9THJOy3czxx5fwM336yJuquGhR8v/+X5kC/pvfJI94RJkyPnFiOdb3eNKkZmYN3H57aQT2ox+V8/z55frUqatD4NOfbr039CDBD+iMZcvKupRjjy1rVXiwX/4yedKTyrSpV7xiw8+99dbSPfWoo0ydheFYurSsifvKV8rH06eXoHfccaWRSr+o694KrnVdppRefHE5Lr20BNWxY5MnPKGEwCc8oYTCPfc0NRRaNtzg18NzIYCeMG5cGem74ILSoVIXuAebObP8TI46auPP3X335M1vLm3fL7usBEZg3e65pwS8Sy9NzjijTKnec8+2q2pGL4W+pNRz8MHleOMbSwC/7LLVQfBd71rdMGe77ZIpU0oIXPswOgg9xYgfsHEzZ5YpVT/+cVn7QVHXZQRvv/1K44bhWLiwfM2ee5bRwl57wwe94E9/Kuve/vjHMtNgY6PpdNfdd5dmVbNnl+O668r5j38s09+H7LprCYBPeUryD/+wZesRgfUy1RPonEWLyoL/E08sXekofv/7ckf8c58rTXCG68tfLj9L22TAQ/3qV6VT5bJlpYnL057WdkUM17JlpXnNUCCcPTuZNSu5/PLyN+S9701e85rhNRYDhk3wAzrrRS8qf7xvvtl6jiHve1/yjneUtXu77Tb8r1u1qnT4vPvu8qZIt1QoLrggefnLSzOT73+/v9bxjWa//nWZMvqznyXTpiUf/7jZI9BBww1+3r0Bw3P88cm8eas76lGmwD7hCZsW+pISnD/60WTu3OSTn2ymNhhpPv3p0kTqoIPKTSahr3889rFl641vfSt54IHkWc9KXvjC5Prr264MRhXBDxieo44qTUzOP7/tSnrD3LnlLvaxx27e1x9xRGma8/73r26bDqPRqlXJm96UvP71yTHHJD/9aekkTH+pqjJz5A9/SD74wfK/80EHlb1N//zntqvbfCtWlLXbMAKY6gkM35FHJjfcUO7SjvamJJ/6VPL3f1+aGkyZsnmvMXt2WSO4alWy//7lTdCax9SpyTbbdLbuXrV4cfKLX5RmEAcf3HY1dMuiRcmrXlVuKL3+9cnHPmb912gxf37y9rcnZ52VPPzhZf3fSSeVLSNGgnvuSb7whTJSPW9eWcO4zz7Jvvs+9Lz33qb00yhr/IDO+7d/S1772rLZ8CGHtF1Nu57xjGTBguSaa7bsdS67LPnhD0ujmN//vgTrlSvL57baKpk8OTnwwNVhcNKk8sZo7Njy+aHH6/p4/PhkwoQt/7c2oa7LDYQf/KD8+y+5JFmypHzu+OOTd79bAOx3d9xRRvj+539K4HvDG9quiDZcfXVZ/3fJJeV33NveVjazH46xY8sNst12697NyD/+MfnEJ5Kzzy4jfc96Vvl7MHducuONpSPtjTeWRjdr2nXXEgT33rv8Hp80qdS95uMJE9xUZbMIfkDn3X57+eP0zneWY7S6884yFe2tby13qTtp6dIyijgUBNcMhGu2SR+ugYHypmL33dd/7Lprd+6yP/BA8pOflKD3wx+WN0hJeeN25JHJc56T/Pd/l8YPDzyQvPSlZb+wAw5ovja6Z8WK0rnzla8sIyXnnFP262P0quvkO98pWz7MmbPpX7/zzuVm5KMfvfp41KM6O2PissvK2uxvf7vcZHvFK0pgffSjH/rcVavKiOZQCFzzfNNNyW23ldHutY0f/9BAuMceyV57lS2A9tyz/O3RYI21CH5AM5761OTee8uo32g1tB3DFVeU7pzdsHRpmRq6YEEZEVyxohwbenz//aXj6JrHvHnJ8uUPfu0xY5KJE0tIHD/+wcf22z/02vjxybhxJSxuvfX6j7Fjy93ryy8vQe/nPy/fe/z40tHvyCOT5z637IO4prvvTj7ykTKddvHi5IQTSvfUyZO787Oms+64o+xZefnl5fyrX5U3vRMnJhdemDz+8W1XSK9YujS56qrVsx42ZsmSsmbwN78pxzXXrJ45MHZsuWk0FAQPOaRMu5w0qfyuG44VK0rQ++hHy/9/H/awsnXP615XXmdz1XX5/TxvXjluu+2hj2+7rfzOXrz4wV+79dblht1QENxzzxIM99uvjD6OlKmydJTgBzTj4x8vjRiuv370vhF/4QvL9KQbbxx503JWrSojlrfc8uBAePvtZZRt4cJyLFq0+vGax+aMOiZl2ubznlfC3pOfPLw78QsWJB/6UNk7ctmy5K//ukwD23ffzauB5i1fXt6Arxn0hkZ2x45NHvOY5IlPLN1wn/WsEv6gU1asKLMjhoLg0HHrrQ9+3g47rB5RG5oRsebHu+ySfO97pevyn/5U1mC/8Y3Jq19dblx1S12XG2E331yOuXNXPx46brml/LuT5G/+JvniF7tX32iyfHnZjmTOnOTkk9uu5iEEP6AZN95Y3nh/6EPJP/5j29V038KFZVrRySePvq0Y6rrckV+0qASx5cvLG47ly9d/rFhRQt8ee2z+950/v3QB/Pzny0jA3/xN8s//XO500xseeKCMgn/3u6tHXCZNWh3ynvjE0tJfgwvacNddZTRw7tzVI2vz5q2eBbGumRBJ8pd/WW50HnNM7zYdWrWq3Lj72MfKTIl///fkNa9pu6r+sGRJ8uMfJ+edV/YYvfvusv70ttvKyGsPEfyA5hx2WJnq98tftl1J9513XvLiF5dW5E9/etvVjC633pp84AOlyVBVlTWBW29d3vjU9erzuh5vvXUJn3vvXaZF7bVXebz77j33B3zEuf/+5PnPL2ugTj01ecpTStDbY4+RNyLO6FTXJRwOhcHbbis3rA4/vO3Khm/lyjKr4mc/Kx2Su7UMod/cf39pOnb++WXU94EHkh13LNsvHX98WZ6w/fZtV/kQgh/QnH/5lzLl7tZbt2ydw0j0yleW9Wrz51tL0Za5c8v+h7/4RQkWVVXWKW7o8dKlZVrUggUPfq2qKv8fXjMQ7r9/MmNG+WPPht17b5m++6tfJeeem7zkJW1XBKPXnXeWwFdVyZVXDr876mh3991lvfF55yU/+lH5e7HLLmWf3uOPL11bx41ru8oNEvyA5lx7bdli4LOfLXf4R4vly8uapOOOS770pbarYXMsXrx6rcxNN5Xzmo9vvrlMY33EI0rX1lNPTbbdtu2qe9M995S731ddlXz96+UNEtCuX/2qTFF9xjPKiFWvTlHthrpO/vznMhV2/vzV57Uf//a3ZcR0zz3L77Hjjy9r0UfQz07wA5p1wAFlKtePf9x2Jd1z8cVleuF3vlPWfNB/Vq1Kfv3rsobwRz8qI4DveU8Z6R1BbwIad9dd5b+F3/0u+da3/PcAvWRoz923v738/hotliwpv4/+4z/Ktkjz5z90P8WkzNZ55CPLseuupePri160erR0BBpu8DNPCdg8xx9fGrzcddfomU4yc2aZ2//sZ7ddCU0ZMyaZPj256KLkP/8z+ad/Kp38PvKRMr306KNH7BuDjrnzztKRc9as0ur++c9vuyJgTSedVLrqvve9ZZ3i0Ue3XVGzZs1KvvCF5CtfKdM299+/rMHfddfV4W7N88MfPmr3QjTiB2yeK65IHve4MuXx1a9uu5rmrVpVRjif+MSyDoDRoa7LHeR//ueyhcmTn5z867+W82h0xx0l9F1/fRn5fs5z2q4IWJfFi8vvqT/9qfy93n//tivqrCVLSgOWL3yhNLQZO7Ysw3jta8s011EW7IY74je6fipA5xx2WJkPP3Nm25V0x69+VTq9HXdc25XQTVVVGpb8/vdlO4k5c8r6mRe+sFwbTebPL2+obrihbNsg9EHv2m67cpOyqso0xkWL2q6oM2bPTk4/vdyIPeGEso/hBz9Yzt/4RvLMZ4660Lcp/GSAzVNVJQRddFFpd9zvZs4sdxSPOqrtSmjD1luXO8k33FC62l5ySVkXcuKJyTnnlPWAv/51aRDTL2+w1jRvXpk6ddNNpdX5M5/ZdkXAxuy7b/K1r5XmJX/7t2UGw0hT16Xp1le/Wm48HXBA8qlPlccXX1xmH7zlLWUKJxtlqiew+S69tLwZ/MY3+ruNe12XPzZ7713e4MNdd5U9BT/zmdL6e23bbpvsvPNDj2nTyvYHe+3V/Zo31803J0ccUUb8fvCDMuIJjBzvfnfyrnclZ56ZnHJK29Vs2O23lxk2V1yx+nzHHeVz++5b1i+eeGJZr8f/0tUTaN7Klcluu5W7/+ee23Y1zfnDH5KDDhp921ewcfffX0bD7rzzocdddz344zvuKPveJcmjHlUC4JFHJk996qZvGXHzzWUfw6Hj1ltLqDzssOSxjy3n/fbb8ilPN91U7qzfdVfZv/KJT9yy1wO6b9WqsgH5xRcnP/958vjHt11RcffdZb/BNYPeLbeUz40ZU35PTp9e+gk87nHlsWmc6yT4Ad1x0kllD68FC5Jttmm7mma87nWlPfZNN5WgC5ujrkv3uR/+sByXXlpGC7fbroSroSA4efKDO4euWFGmag2FvMsuK8EvScaPL2/i9tqrbK3w29+W/SaTsgH9Yx6zOggedlgyZcrqN051XaZp33bb6mP+/Ad/fPXV5ftfdFHpDgiMTHffXYLT8uUlbO2ySzt1rFpVAugnP1lmEAyZMmV1uJs+vfzuGhhop8YRSPADuuMHPyjt3L/3vf5s63777ck++5RF5F/8YtvV0E8WLizhbygIXn99ub7ffiUAPuIRJeRdfnl5blIaGjz5ycmTnlTOj350WXs6ZNmy0nTmyivLmsMrr0x+85vV01EHBpKpU8umxrfdtu71iOPGlRscu+1Wvt9b31rehAEj21VXld8dT3xiWbYwtou7uj3wQNlf79OfLg1adt01ec1rynKRww5Ldtqpe7X0IcEP6I6lS5OJE5OXvrQ/g9EZZ5T9CmfNKnckoSlz5pSRtR/8IPnJT0q78kMOKQFv6NictYHLlyfXXrs6CF53XdnHaijcDR277lrOD3uYvQqhX335y2WN3PTpZT/eo49ODj64uf/m//jHskzirLPKVPfHPS75+78vfQHGjWvme45Cgh/QPa94RZm6MXdumbbWL/7859LQ5cgjy3RW6JalS0tgM9UJ6LQzz0zOPrvcCErKDaWjjy7HM56x6WuO11bXyU9/WrpvXnBBstVWyYtfnLz+9ckTnuDGUgMEP6B7vvOd5Nhjy927ww5bPQ3tSU8a2S2W3//+snH3VVclhx7adjUA0Dnz5iXf/37Zl/Pii8vU7+23T571rBICjzoqmTRpw6+xcmWZxnn//cl995V1yJ/6VHLNNWW6+mtfW5qi7b57d/5No5TgB3TXD39YpqdddlnpzLVsWbm+//4PXpN04IEjoyvXokVlbd/06eUPIwD0qyVLyv6kF15YguDcueX6Yx9bumsOBbs1z/ffv3r98ZoOOaRM53z5y/trFlAPE/yA9ixdWtYTrdmFcGgfnh13LAHw058uTSx61ac/Xaal/Pzn9i0DYPSo6zJi993vlmP+/GTChGSHHcp5Q4/32ad0Gjads6sEP6B31HVpXHHZZSUIfvWryXHHlXMvWrastNTfe+8S/AAAetRwg18X+7gCo1ZVlSA1eXLyV39V9vv7wheSj32svb2ENuScc8o+aV/4QtuVAAB0xAhYaAP0nVNOKaNqZ53VdiUPtXJl8sEPlmYuRx7ZdjUAAB0h+AHd96hHJUcckXz+8yVo9ZJvf7tsLnvGGdYoAAB9Q/AD2nHqqaVrWC91zKzr5AMfKBu1v+hFbVcDANAxgh/QjmOOKfsDnXlm25WsdvHFZUPbN7+5bDgLANAnBD+gHVtvnZx8ctn/b86ctqspPvCBssnsq17VdiUAAB0l+AHtOemkMrL2+c+3XUnyy1+WzWtPP710HQUA6COCH9CeSZPKfn5nn50sXtxuLR/4QPKIR5QwCgDQZwQ/4P+3d/+xd9X1HcefL9oijE6wUjYilIoQpwmlSld1MmBAtELVaRzSDIONgtIaRYcLEKPBRCNzgRlFIqMFYlhBEBSRwNS5gFnGWrCWMmagKtjUUaUTbBkI63t/nEOsTX/zvffcc3k+kpt7zud7ez7v5J17vn197/mc263Fi2HDBrj++u5quO8++Na34MMfhqlTu6tDkiRpQAx+krp1/PHN1zt0eZOXz32uCXwf+lB3NUiSJA2QwU9St5Lmqx2WL28ew7ZmDVx3HXzwgzBt2vDnlyRJGgKDn6Tuvec9sN9+cPnlw5/785+HyZPhYx8b/tySJElDYvCT1L3994czzoBly5r1fsOybh1cdRUsXAgHHzy8eSVJkobM4CdpNCxaBE89BVdfPbw5L70Unn0WPv7x4c0pSZLUAYOfpNEwaxYce2xzuefmzYOfb8OGZq7TT4dXvGLw80mSJHXI4CdpdCxaBA89BN/97mDnqYILL4RNm+D88wc7lyRJ0ggw+EkaHe98Jxx0EFx22WDn+exn4StfgfPOg6OOGuxckiRJI8DgJ2l0vOhF8P73w623wsMPD2aOK6+ET3yiuZnMxRcPZg5JkqQRY/CTNFrOPrt5vuKKiT/2LbfABz4A8+bB0qWwl6dASZL0wuD/eiSNlsMOg/nzm0/mnn564o77gx/Au98NxxwDN9wAU6ZM3LElSZJGnMFP0uhZvBjWr4ebbpqY461eDW99K8yYAd/+NkydOjHHlSRJ6gmDn6TRc/LJcMQR8OUvP/9jPfJIcytZUYkAAAi5SURBVGnnvvvCHXfA9OnP/5iSJEk9Y/CTNHr22gvOOae5PHPVqj0/zmOPwZvfDBs3wu23w8yZE1aiJElSnxj8JI2m974X9tmn+ZL1PbFpE5x6Kvz0p81NXWbNmtDyJEmS+sTgJ2k0TZsGCxbAV78Kjz++e//2mWfgtNNg+XJYtgyOO24wNUqSJPWEwU/S6Fq0qPnkbuFCWLKkCXJPPrnjf1MFZ50Ft93WrBF8xzuGU6skSdIIm9x1AZK0XXPmNN/rd+21cPPNzVgCRx7ZXLo5axYcdVTzPHNmszbwggvgmmvgooua7+yTJEkSqaqua9hjc+bMqRUrVnRdhqRB27y5Wau3atXvP9asaT7hg+YrGo44AlaubG4Mc9llTUiUJEkaY0nuqao5O32dwU9Sb23cCPff//thcPZsuOQSmDSp6+okSZIGbleDn5d6SuqvqVPhda9rHpIkSdoub+4iSZIkSWPO4CdJkiRJY87gJ0mSJEljzuAnSZIkSWPO4CdJkiRJY87gJ0mSJEljzuAnSZIkSWPO4CdJkiRJY87gJ0mSJEljzuAnSZIkSWNuYMEvydIk65Os3mJsdpJ/T7IyyYokc9vxE5I83o6vTPLJQdUlSZIkSS80g/zE72pg3lZjfwdcVFWzgU+2+8+5q6pmt49PD7AuSZIkSXpBGVjwq6o7gQ1bDwMvbrf3B9YNan5JkiRJUmPykOc7F7gjyd/ThM4/2+Jnb0jyI5oweF5V3b+tAyQ5GzgbYMaMGQMuV5IkSZL6b9g3dzkH+GhVHQp8FFjSjt8LHFZVRwNfBL6xvQNU1RVVNaeq5kyfPn3gBUuSJElS3w07+J0J3NRu3wDMBaiqJ6pqY7t9GzAlyYFDrk2SJEmSxtKwg9864Ph2+0TgQYAkf5wk7fbctq7HhlybJEmSJI2lga3xS7IMOAE4MMla4FPAWcAXkkwGnqJdqwe8CzgnybPA/wKnV1UNqjZJkiRJeiEZWPCrqgXb+dEx23jtl4AvDaoWSZIkSXohG/alnpIkSZKkITP4SZIkSdKYM/hJkiRJ0pgz+EmSJEnSmEufb56Z5JfAw13XsQ0HAr/qugjtMfvXb/av3+xfv9m/frN//Wb/+u359O+wqpq+sxf1OviNqiQrqmpO13Voz9i/frN//Wb/+s3+9Zv96zf712/D6J+XekqSJEnSmDP4SZIkSdKYM/gNxhVdF6Dnxf71m/3rN/vXb/av3+xfv9m/fht4/1zjJ0mSJEljzk/8JEmSJGnMGfwmWJJ5SX6c5KEk53ddj3YsydIk65Os3mJsWpLvJHmwfX5JlzVq25IcmuT7SR5Icn+Sj7Tj9q8HkuyT5D+S/Kjt30Xt+MuT3N327/oke3ddq7YvyaQkP0xya7tv/3oiyc+S3JdkZZIV7Zjnz55IckCSG5P8V/t78A32rx+SvLJ93z33eCLJucPon8FvAiWZBFwGvAV4NbAgyau7rUo7cTUwb6ux84HvVdWRwPfafY2eZ4G/qapXAa8HFrfvN/vXD08DJ1bV0cBsYF6S1wMXA5e2/fsf4H0d1qid+wjwwBb79q9f/qKqZm9xC3nPn/3xBeD2qvoT4Gia96H964Gq+nH7vpsNHAM8CdzMEPpn8JtYc4GHquonVfVb4Drg7R3XpB2oqjuBDVsNvx24pt2+BvjLoRalXVJVv6iqe9vt39D80nsZ9q8XqrGx3Z3SPgo4EbixHbd/IyzJIcCpwJXtfrB/fef5sweSvBg4DlgCUFW/rapfY//66CRgTVU9zBD6Z/CbWC8Dfr7F/tp2TP3yR1X1C2jCBXBQx/VoJ5LMBF4D3I396432MsGVwHrgO8Aa4NdV9Wz7Es+ho+0fgL8FNrf7L8X+9UkB/5zkniRnt2OeP/vhcOCXwFXtpdZXJtkP+9dHpwPL2u2B98/gN7GyjTFvmyoNUJKpwNeBc6vqia7r0a6rqv9rL3U5hOaKiVdt62XDrUq7Isl8YH1V3bPl8DZeav9G1xur6rU0y1MWJzmu64K0yyYDrwUur6rXAJvwss7eaddAvw24YVhzGvwm1lrg0C32DwHWdVSL9tyjSQ4GaJ/Xd1yPtiPJFJrQd21V3dQO27+eaS9R+leatZoHJJnc/shz6Oh6I/C2JD+jWdZwIs0ngPavJ6pqXfu8nmZ90Vw8f/bFWmBtVd3d7t9IEwTtX7+8Bbi3qh5t9wfeP4PfxFoOHNne1Wxvmo9vb+m4Ju2+W4Az2+0zgW92WIu2o11PtAR4oKou2eJH9q8HkkxPckC7vS9wMs06ze8D72pfZv9GVFVdUFWHVNVMmt91/1JVf43964Uk+yX5w+e2gTcBq/H82QtV9d/Az5O8sh06CfhP7F/fLOB3l3nCEPrnF7hPsCSn0PzVcxKwtKo+03FJ2oEky4ATgAOBR4FPAd8AvgbMAB4B/qqqtr4BjDqW5FjgLuA+frfG6EKadX72b8QlmUWzeH0SzR8hv1ZVn05yOM0nSNOAHwJnVNXT3VWqnUlyAnBeVc23f/3Q9unmdncy8E9V9ZkkL8XzZy8kmU1zY6W9gZ8AC2nPpdi/kZfkD2juC3J4VT3ejg38/WfwkyRJkqQx56WekiRJkjTmDH6SJEmSNOYMfpIkSZI05gx+kiRJkjTmDH6SJEmSNOYMfpIk7aYkG7fYPiXJg0lmdFmTJEk7MrnrAiRJ6qskJwFfBN5UVY90XY8kSdtj8JMkaQ8k+XPgH4FTqmpN1/VIkrQjfoG7JEm7KckzwG+AE6pqVdf1SJK0M67xkyRp9z0D/Bvwvq4LkSRpVxj8JEnafZuB04A/TXJh18VIkrQzrvGTJGkPVNWTSeYDdyV5tKqWdF2TJEnbY/CTJGkPVdWGJPOAO5P8qqq+2XVNkiRtizd3kSRJkqQx5xo/SZIkSRpzBj9JkiRJGnMGP0mSJEkacwY/SZIkSRpzBj9JkiRJGnMGP0mSJEkacwY/SZIkSRpzBj9JkiRJGnP/D1zMR8L3fpUOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x720 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(15,10))\n",
    "plt.xlabel(\"K\")\n",
    "plt.ylabel(\"Error rate\")\n",
    "plt.plot(error_rate, \"r\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "20"
      ]
     },
     "execution_count": 114,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.argmin(np.array(error_rate))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SVD-based CF MSE: 185.15160022347177\n"
     ]
    }
   ],
   "source": [
    "#get SVD components from train matrix. Choose k.\n",
    "u, s, vt = svds(train_data_matrix, k = 20)\n",
    "s_diag_matrix=np.diag(s)\n",
    "X_pred = np.dot(np.dot(u, s_diag_matrix), vt)\n",
    "print('SVD-based CF MSE: ' + str(mse(X_pred, test_data_matrix)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.42253129e-01,  1.91892606e-01,  4.07786483e-01, ...,\n",
       "         4.78271307e-02,  1.83263808e-02,  2.02365759e-02],\n",
       "       [ 1.37099162e-02,  2.94124351e-02,  5.66586449e-02, ...,\n",
       "        -1.63267887e-03,  1.85948728e-03,  6.53201990e-03],\n",
       "       [ 4.85809658e-03,  1.32983418e-02,  2.44914582e-02, ...,\n",
       "         2.21900761e-03,  6.02568118e-04,  4.43749642e-03],\n",
       "       ...,\n",
       "       [ 6.47250542e-03,  1.56418829e-02,  3.04814204e-02, ...,\n",
       "        -2.04059891e-03,  8.46335469e-04,  3.97239789e-03],\n",
       "       [ 3.03265186e-02,  6.96778113e-02,  1.25011698e-01, ...,\n",
       "         7.69204833e-03,  2.64959357e-03,  2.67728361e-02],\n",
       "       [ 3.26790011e-01,  7.04080096e-01,  1.29399046e+00, ...,\n",
       "         7.69769752e-03,  5.12890987e-02,  1.77526430e-01]])"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[12.,  0.,  1., ...,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0., ...,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0., ...,  0.,  0.,  0.],\n",
       "       ...,\n",
       "       [ 0.,  0.,  0., ...,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0., ...,  0.,  0.,  0.],\n",
       "       [ 3.,  3.,  0., ...,  0.,  0.,  0.]])"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_data_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 推荐前20个"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_top_items(items_dictionary, values):\n",
    "    ids = []\n",
    "    for i in range(0, len(values-1)):\n",
    "        ids.append(list(items_dictionary.keys())[list(items_dictionary.values()).index(values[i])])\n",
    "    \n",
    "    return ids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "ids = get_top_items(items_dictionary, np.argsort(user_prediction_non_biased[3])[::-1][:20])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "ids_items = get_top_items(items_dictionary, np.argsort(item_prediction_non_biased[3])[::-1][:20])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "ids_svd = get_top_items(items_dictionary, np.argsort(X_pred[3])[::-1][:20])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## user_top_recmandation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'5a905f000fc1ff3df7ca807d57edb608863db05d'"
      ]
     },
     "execution_count": 123,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    " list (users_dictionary.keys()) [list (users_dictionary.values()).index (user)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "总共生成 98632 次推荐。\n"
     ]
    }
   ],
   "source": [
    "import pickle\n",
    "j=0\n",
    "recomList={}\n",
    "for user in users_dictionary.values():\n",
    "    pred = get_top_items(items_dictionary,np.argsort(user_prediction_non_biased[user])[::-1][:20])\n",
    "    user_id=list (users_dictionary.keys()) [list (users_dictionary.values()).index (user)]\n",
    "    recomList[user_id]=pred\n",
    "    \n",
    "    #print(' user: %s,  song: %s' % (user_id,pred))\n",
    "    j+=1\n",
    "print (\"总共生成\",j,\"次推荐。\")\n",
    "output = open('recommendationList_user.pkl', 'wb')\n",
    "pickle.dump(recomList, output)\n",
    "output.close()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Item top recommandation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "总共生成 98632 次推荐。\n"
     ]
    }
   ],
   "source": [
    "import pickle\n",
    "j=0\n",
    "recomList={}\n",
    "for user in users_dictionary.values():\n",
    "    pred = get_top_items(items_dictionary,np.argsort(item_prediction_non_biased[user])[::-1][:20])\n",
    "    user_id=list (users_dictionary.keys()) [list (users_dictionary.values()).index (user)]\n",
    "    recomList[user_id]=pred\n",
    "    \n",
    "    #print(' user: %s,  song: %s' % (user,pred))\n",
    "    j+=1\n",
    "print (\"总共生成\",j,\"次推荐。\")\n",
    "output = open('recommendationList_item.pkl', 'wb')\n",
    "pickle.dump(recomList, output)\n",
    "output.close()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## svd top recmmandation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "总共生成 98632 次推荐。\n"
     ]
    }
   ],
   "source": [
    "import pickle\n",
    "j=0\n",
    "recomList={}\n",
    "for user in users_dictionary.values():\n",
    "    pred = get_top_items(items_dictionary,np.argsort(X_pred[user])[::-1][:20])\n",
    "    user_id=list (users_dictionary.keys()) [list (users_dictionary.values()).index (user)]\n",
    "    recomList[user_id]=pred\n",
    "    \n",
    "    #print(' user: %s,  song: %s' % (user,pred))\n",
    "    j+=1\n",
    "print (\"总共生成\",j,\"次推荐。\")\n",
    "output = open('recommendationList_svd.pkl', 'wb')\n",
    "pickle.dump(recomList, output)\n",
    "output.close()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## precision&recall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pprint, pickle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [],
   "source": [
    "#测试集转为字典\n",
    "user_song = dict()\n",
    "for _, row in test_data.iloc[:,0:3].iterrows():\n",
    "    user, song, listen_count = row\n",
    "    user_song.setdefault(user, {}).update({str(song): listen_count})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getRecom(rec_list='recommendationList_user.pkl'):\n",
    "    pkl_file = open(rec_list, 'rb')\n",
    "    recomList = pickle.load(pkl_file)\n",
    "    #pprint.pprint(recomList)\n",
    "    pkl_file.close()\n",
    "    return recomList"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [],
   "source": [
    "#precision-前20个推荐物品中用户喜欢物品的比例\n",
    "def calPrecision(recom,user_item):\n",
    "    precisionForUser={}  #precision for every user\n",
    "    for key in recom:\n",
    "        hit_Num=0\n",
    "        if key in user_item:\n",
    "            for i in range(20):\n",
    "                item=recom[key][i]\n",
    "                if item in user_item[key]:\n",
    "                    #如果用户u的推荐列表中每有一个item 在 test中用户u的记录中出现，hit_num+1\n",
    "                    hit_Num+=1\n",
    "            precisionForUser[key]=hit_Num*1.0/10\n",
    "    Precision=sum(precisionForUser.values())*1.0/(len(precisionForUser))\n",
    "    return Precision"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [],
   "source": [
    "#recall-推荐命中的物品数量：用户看过的数量\n",
    "def calRecall(recom,user_item):\n",
    "    recallForUser={}          #recall for every user\n",
    "    for key in recom:\n",
    "        hit_Num=0\n",
    "        if key in user_item:\n",
    "            for i in range(20):\n",
    "                item=recom[key][i]\n",
    "                if item in user_item[key]:\n",
    "                    #如果用户u的推荐列表中每有一个item 在 test中用户u的记录中出现，hit_num+1\n",
    "                    hit_Num+=1\n",
    "            recallForUser[key]=hit_Num*1.0/(len(user_item[key]))\n",
    "    Recall=sum(recallForUser.values())*1.0/(len(recallForUser))     \n",
    "    return Recall"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## user based CF-evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [],
   "source": [
    "recom=getRecom(rec_list='recommendationList_user.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [],
   "source": [
    "precision=calPrecision(recom,user_song)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [],
   "source": [
    "recall=calRecall(recom,user_song)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 该推荐算法准确率为: 0.10049967320953868,  召回率为: 0.0893159581529443\n"
     ]
    }
   ],
   "source": [
    "print(' 该推荐算法准确率为: %s,  召回率为: %s' % (precision,recall))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## item based CF-evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [],
   "source": [
    "recom=getRecom(rec_list='recommendationList_item.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [],
   "source": [
    "precision=calPrecision(recom,user_song)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [],
   "source": [
    "recall=calRecall(recom,user_song)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 该推荐算法准确率为: 0.004345259429487271,  召回率为: 0.0070275225038311795\n"
     ]
    }
   ],
   "source": [
    "print(' 该推荐算法准确率为: %s,  召回率为: %s' % (precision,recall))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## svd based CF-evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [],
   "source": [
    "recom=getRecom(rec_list='recommendationList_svd.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [],
   "source": [
    "precision=calPrecision(recom,user_song)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [],
   "source": [
    "recall=calRecall(recom,user_song)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 该推荐算法准确率为: 0.09456262781728206,  召回率为: 0.09013356561323906\n"
     ]
    }
   ],
   "source": [
    "print(' 该推荐算法准确率为: %s,  召回率为: %s' % (precision,recall))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
